Method for identifying a nuclear explosion based on krypton and xenon isotopes

ABSTRACT

The invention pertains to the field of nuclear physics and can be used in system for identifying nuclear explosions based on the measured activities in the atmosphere of naturally-occurring radioactive gases (NORG). The technical result is an increase in the determination efficiency and in the reliability of punctual estimations of deposits from various types of fission in the global activity for each krypton and xenon isotope.

FIELD OF THE INVENTION

The invention relates to nuclear physics and may be used in systems for identifying radioactivity in the atmosphere.

PRIOR ART

The identification of nuclear explosions based on isotopes of radioactive noble gases (RNG) takes place in the process of monitoring radioactive conditions for the purpose of monitoring the compliance with the Comprehensive Nuclear-Test-Ban Treaty (CTBT).

A method for distantly detecting nuclear material is known, which comprises irradiating an intercepted object with an electron beam and registering neutrons radiated from the object (U.S. Pat. No. 4,320,298, 1982).

Disadvantages of the said method are its non-applicability in the atmospheric conditions and its non-applicability for detecting objects that do not contain deuterium.

Methods for distantly detecting nuclear materials are known, which comprise irradiating an intercepted object with a neutron beam and registering neutrons or gamma-rays radiated from the object (SU Inventor's Certificate # 439740; SU Inventor's Certificate # 1349478; U.S. Pat. No. 4,483,817).

A disadvantage of these methods are their non-applicability in the conditions where for some or other reasons external radiation effect on an intercepted object is prohibitive.

A method for distantly detecting nuclear materials is known, which comprises detecting gamma rays intensity in the energy range from 0.1 to 2.0 MeV near an intercepted object (Sagdeyev, R. Z. et al. Problems of Monitoring sea-based cruise missiles with nuclear warheads. Preprinted by IKI of the USSR AS, Pr-1373.-M., 1988.).

A disadvantage of this method is the possibility of false detection of nuclear material, since radiation with such energy may be also created by non-explosive objects containing radioactive materials.

Methods for identifying nuclear explosions by radioactive isotopes of krypton and xenon are known.

A method for identifying nuclear explosions by radioactive isotopes of krypton and xenon (SU Inventor's Certificate # 366771) is known. The mode of measuring RNG activity in this method is as follows: a sample is taken in the atmosphere (after an event), which is studied for some time. During this time the krypton and xenon isotopes activity is measured with a single-crystal scintillation gamma spectrometer NaJ(Tl).

Thus measured activities of krypton and xenon isotopes are used for forming a system of linear algebraic equations (SLAB) relative to unknown contributions of RNG sources to total activity of krypton and xenon isotopes.

The equation system is solved with the use of least-squares method (LSM).

The closest as to the technical essence is a method that enables to eliminate some of the disadvantages of the above method, i.e., to take into account the matrix element errors in an equation system to be solved and ensure obtaining of a stable solution with the use of the regularization method by A. N. Tikhonov (Greshilov A. A., Tetjukhin A. A. “An Algorithm for Identifying Sources of Radioactive Noble Gases”. Bulletin of MGTU Named After N. E. Bauman. Natural Sciences Series, 2003. No. 2, p. 3-19.)

The known method for identifying nuclear explosions by radioactive krypton and xenon isotopes comprises:

1. Measuring activities Ã_(i)(t) of separate krypton and xenon isotopes in the atmosphere (t is measurement time), where i= 1,n, n is a number of measured isotopes.

2. Determining a time interval [t_(H),t_(k)] in which separation occurs, for different types of fissionable materials (a fission type is understood as one of the fission variants of heavy nuclei of U²³⁵, U²³⁸, PU²³⁹ by neutrons of the fission spectrum or neutrons with energy of 14 MeV) over activity relations of krypton and xenon isotopes, which are built with and without due regard to separation;

3. Setting time grid with At spacing on the [t_(H),t_(K)] interval;

4. Forming and memorizing two-dimensional signal {a_(ij)(t_(q),t)}, i= 1,n, j= 1,m of specific activities of krypton and xenon isotopes for each grid node t_(q) in which lines correspond to particular isotopes and columns to a particular fission type; n is a number of isotopes studied; m is a number of fission types studied; t is time of sample measuring; t_(q) is supposed separation time;

5. Regarding measured activity values A, and elements of the two-dimensional signal {a_(ij)(t_(q),t)} as statistically independent quantities distributed according to the normal law with root-mean-square deviations σ(A_(i)(t)) and σ(a_(ij)(t_(q),t)), respectively, a regular signal is formed at the set time t_(q):

${F_{l} = {\frac{1}{2}{\sum\limits_{i = 1}^{n}\begin{pmatrix} {\frac{\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)} \cdot \left( {\rho \; N_{j}} \right)}}} \right)^{2}}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)} +} \\ {\sum\limits_{j = 1}^{m}\frac{\left( {{a_{ij}\left( {t_{q},t} \right)} - {a_{ij}^{true}\left( {t_{q},t} \right)}} \right)^{2}}{\sigma^{2}\left( {a_{ij}\left( {t_{q},t} \right)} \right)}} \end{pmatrix}}}},$

where:

(ρN_(j)), j= 1,m

a_(ij) ^(true)(t_(q),t) are contributions of radioactivity sources to a total activity;

a_(ij)(t_(q),t) are unknown true values of specific activity;

a_(ij)

(t_(q),t) are specific activities calculated over independent and cumulative elements, which have errors, of radioactive transformation isobaric chains (RTIC);

l is iteration number for search of estimates for {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} and â_(ij) ^(true)(t_(q),t).

6. Setting the numbers of γ1, γ2 that characterize the accuracy of estimates of {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} and â_(ij) ^(true)(t_(q),t).

7. Finding iteratively a minimum signal F1 from {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} and â_(ij) ^(true)(t_(q),t) by using the regularization by A. N. Tikhonov and complying with the limitations |a_(ij)(t_(q),t)−â_(ij) ^(true)(t_(q),t)|≧3σ(a_(ij)(t_(q),t)), i=1, 2, . . . , n, j=1, 2, . . . m , until the following conditions are met:

${{\frac{{()}_{l} - {()}_{l - 1}}{{()}_{l}}} < \gamma_{1}},{{\frac{F_{l} - F_{l - 1}}{F_{l}}} < \gamma_{2}},$

where the index (l−1) is for the value obtained in the preceding iteration.

8. Estimate covariance matrix is determined by the relation

${D_{ij} = \left( {- \frac{\partial^{2}F_{l}}{{\partial{()}}{\partial{()}}}} \right)^{- 1}},i,{j = 1},2,\ldots \mspace{14mu},{m.}$

But, this known method has the following disadvantages: 1) only the Tikhonov's regularization is used, which requires the regularization parameter to be determined additionally, and there is no unique method for finding it; the Tikhonov's regularization “smoothes” a solution, which may cause great errors in identification;

2) the considered equation system is overdefined, i.e., the number of lines of a two-dimensional signal {a_(ij)(t_(q),t)} is greater than that of columns−the number of equations is greater than that of unknowns. A possibility of applying this approach in real situations, when radioactivity of krypton and xenon isotopes is measured in several days after the event, and when the number of unknown ρN_(j) (contributions from various sources) is greater than the number of isotopes to be measured, i.e., the number of Ã_(i)(t) quantities, has not been studied.

Thus, the known solution does not ensure identification of a nuclear explosion in the most probable case, when activity of 2-4 isotopes is measured, and samples are taken in 5-6 days after a corresponding event.

The said disadvantages, apparently, will not enable to apply this method in real conditions due to its low practical efficiency.

SUMMARY OF THE CLAIMED INVENTION

The technical effect of the claimed method is improvement in reliable etermination of a fact of a nuclear explosion made, the number of measured isotopes being less than a considered number of unknowns (fission types).

Efficiency of the claimed method is ensured due to:

1) simultaneous consideration of various supposed combinations of krypton and xenon isotope activity sources and various mechanisms of a nuclear explosion;

2) development and inclusion into the nuclear explosion identification method of a multi-criterion mathematical programming mechanism enabling to take into account all possible types of additional conditions (non-negativity of a solution, boundedness of solution) that a solution estimate should comply with, and elimination of the regularization parameter (by A. N. Tiknonov), clearly formalized procedures that have no definition;

3) combination of two types of U²³⁵ fission by fission spectrum neutrons and by neutrons with energy of 14 MeV in a single fission type as well as two types of Pu²³⁹ fission by fission spectrum neutrons and by neutrons with energy of 14 MeV in a single fission type by averaging independent and cumulative outputs of the RTIC elements.

The technical effect of the claimed invention is achieved by developing a method for identifying a nuclear explosion by radioactive isotopes of krypton and xenon, which is characterized by measuring, at the time t after an event, signals A_(i)(t) describing a change in total activity of each isotope in the atmosphere near a measuring station, building relations between isotope activities and time without due regard to separation and relations between isotope activities which are drawn from measured points in reverse time for all considered fission types, determining a separation interval [t_(H), t_(k)], setting a time grid within that separation interval [t_(H), t_(d)], forming combinations of fission types, calculating, for each grid node t_(q) within the separation interval and each combination, a two-dimensional signal A describing values of “specific” activities of each isotope, depending on separation time t_(q), and measurement time t, and a potential source (of a fissionable material and neutron energy), setting mean square values σ(A_(i)(t)) for errors of measured signals A_(i)(t) and errors σ(a_(ii)(t_(q),t)) of elements of a two-dimensional signal {a_(i;j)(t_(q),t}, setting quantities γ₁,γ₂ defining accuracy of calculations for signal estimations pN_(j) and J₂ signal, identifying a nuclear explosion with the use of formation of a signal J₁ determining accuracy of a solution obtained from the sum of squares of difference between the signals A, (t) and products of lines of a two-dimensional signal {a_(i;j)(t_(q),t} multiplied by values of signals pN_(j) and the signal J₂ determining a type of the pN_(j) signal, forming limitation signals and a target signal (target function) from the signals J₁ and J₂ in the set combinations, and finding, with the use of iteration process and with refining during each iteration, values of the {a_(i,j)(t_(q),t} two-dimensional signal of the elements of the pN_(j) signal, determination of point estimates of activity contributions from each fission type possible to total activity, selecting an optimal combination of fission types by a sum of discrepancy squares.

The method is also characterized by combination of two types of U²³⁵ and Pu²³⁹ fission by neutrons from different energy groups in a single fission type η^(U) ²³⁵

η^(Pu) ²³⁹ , respectively, by summing with weights of independent and cumulative outputs of isobaric chain elements corresponding to different fission types, forming a two-dimensional grid according to weights c₁ and c₂, calculating specific activity for each pair of weight values c1, c2 of the two-dimensional signal elements {a_(i;j)(t_(q),t}_(cc), calculating, according to Item 1, signal estimates (pNj)_(cc), selecting an estimate (pNj)_(cc) for which the value of

$\sum\limits_{i = 1}^{n}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{\left( {a_{ij}^{true}\left( {t_{q},t} \right)} \right)_{c_{1}c_{2}}{()}_{c_{1}c_{2}}}}} \right)^{2}$

is minimal.

Kr^(83m), Kr^(85m), Kr⁸⁵, Kr⁸⁸, Xe^(131m), Xe^(133m), Xe¹³³, Xe¹³⁵ are takes as individual krypton and xenon isotopes.

In order to establish a fact of nuclear explosion by a small number of measured isotopes, two types of U²³⁵ fission (fission by fission spectrum neutrons and neutrons with energy of 14 MeV) are combined in a single fission type, and two types of Pu²³⁹ fission (fission by fission spectrum neutrons and neutrons with energy of 14 MeV) are combined in a single fission type by means of averaging independent and cumulative outputs, as corresponding to the said fission types, of RTIC elements.

Activity contributions from individual sources {circumflex over (ρ)}{circumflex over (n)}{circumflex over (ρn_(j))} are evaluated by means of forming several target functions and using methods of multi-criterion mathematical programming, reducing a multi-criterion task to a single-criterion task with limitations, obtaining a solution to the said single-criterion task with limitations by iteration calculation procedures, determining a contribution of each source to total activity according to {circumflex over (ρ)}{circumflex over (ρ_(j))} values (number of fissions of j-type), i.e., identifying parameters of a nuclear explosion. A combination of RNG isotope sources, which gives a minimum value of the F_(l) signal, is taken as true in the solution process.

Best embodiments of the method for identifying a nuclear explosion by radioactive isotopes of krypton and xenon.

The essence of the claimed invention is explained hereinbelow by means of its description and graphic materials, where:

FIG. 1 shows a change in relative activity of isotopes A(Xe133m)/A(Xe133) for a case of fission of U^(f) ₂₃₅ and Pu^(f) ₂₃₉without separation (solid line and dashed line, respectively) and with due regard to separation from preceding isotopes (line with markers), point l is a relation of measured isotope activities at the time t=12 hours.

FIG. 2 shows a general block diagram of the algorithm used for obtaining estimates for separation time {circumflex over (t)}_(q) and a solution {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))}.

FIG. 3 shows a block diagram for finding a solution to poorly conditioned system of linear algebraic equations with the use of confluent analysis (corresponds to Block 1 of the general block diagram shown in FIG. 2).

FIG. 4 shows a plot of dependence of cumulative outputs of isotopes Xe¹³³ and Xe¹³⁵ on relative contribution of fission spectrum neutrons and neutrons with energy of 14 MeV.

For carrying out the method for identifying a nuclear explosion by radioactive isotopes of krypton and xenon it is necessary to determine the following parameters of nuclear explosion sources on the basis of RNG in the atmosphere:

-   -   contributions of each fission type (fissionable material and         neutron energy {circumflex over (ρ)}{circumflex over         (N)}{circumflex over (ρN_(j))}, j= 1,m) to the total RNG         activity measured in the atmosphere.

As a practically justified assumption for the claimed method, registered signals (of isotope activity) are considered as deterministic, subjected to additive interference, which parameter estimates should be defined.

At instant fission i-number isotope appears as a result of different fission types, and its measured activity Ã_(i)(t) may be expressed as follows [3]:

$\begin{matrix} {{{\sum\limits_{j = 1}^{m}{{a_{ij}\left( {\theta,\eta,\lambda,t,t_{q}} \right)}\rho \; N_{j}}} = {{\overset{\sim}{A}}_{i}(t)}},{i = 1},2,\ldots \mspace{14mu},n,} & (1) \end{matrix}$

where: a_(ij)(θ,η,λ,t,t_(q)) is activity of i-isotope at j-type of fission for one act of decay, which is calculated subject to separation at the time t>t_(q), i.e., specific activity;

θ is vector of parameters characterizing separation of measured isotopes from preceding thereto;

η is vector of independent outputs of isotopes (at j-type of fission);

λ is vector of decay constants;

t is observation time;

t_(q) is supposed separation time of separation of krypton and xenon isotopes from preceding isotopes along chains of radioactive transformations;

ρ is proportion of i-isotope in a sample (ρ value is usually unknown);

N_(j) is number of fissions of j-type.

Until the time of separation t_(q), specific activity a_(nj) ^(is q)(θ, η, λ, t_(q)) is determined by the formula:

$\begin{matrix} {{{a_{nj}^{t_{q}}\left( {\theta,\eta,\lambda,t_{q}} \right)} = \left\{ {\begin{matrix} {{\eta_{n}\lambda_{n}{\exp \left( {{- \lambda_{n}}t_{q}} \right)}} +} \\ {\sum\limits_{p = 1}^{p_{\max}}{\sum\limits_{i_{p} = 1}^{n_{p} - 1}{\eta_{i_{p}}\lambda_{n}}}} \end{matrix}\left\lbrack {\prod\limits_{r_{p} = i_{p}}^{n_{p} - 1}\; {\gamma_{r_{p}}\lambda_{r_{p}}{\sum\limits_{s_{p} = i_{p}}^{n_{p}}\frac{\exp \left( {{- \lambda_{s_{p}}}t_{q}} \right)}{\prod\limits_{\underset{q_{p} \neq s_{p}}{q_{p} = i_{p}}}^{n_{p}}\; \left( {\lambda_{q_{p}} - \lambda_{s_{p}}} \right)}}}} \right\rbrack} \right\}},} & (2) \end{matrix}$

where:

η_(i) is independent output of i-isotope;

n_(p) is studied isotope number along p-branch;

n is maximum member from {n_(p)};

p_(max) is number of chain branches;

(n_(p)−1) is number of isotopes preceding the studied one along p-branch of decay;

γ_(r) _(p) is proportion of r-member of the chain, which is obtained from (r−1) along the p-chain;

λ_(t) _(p) , λ_(r) _(p) , λ_(s) _(p) , λ_(q) _(p) are decay constants of isotopes having respective numbers i_(p), r_(p), s_(p), q_(p) along the p-branch, wherein i_(p)≦r_(p)≦n_(p)−1; i_(p)≦s_(p)≦n_(p); i_(p)≦q_(p)≦n_(p) and q_(p)≠s_(p);

t_(q) is time when instant separation of studied isotope from preceding ones occurs, after which isotope decal goes exponentially with decay constant λ_(n).

After the time of separation isotopes decay according to their decay constants λ_(i):

a _(ij)(θ, η, λ, t, t _(q))=a _(ij) ^(t) ^(q) (θ, η, λ, t _(q))exp(−λ_(t)(t−t _(q)), i=1, 2, . . . , n, j=1, 2, . . . , m,   (3)

where a_(ij) ^(t) ^(q) (θ, η, λ, _(q)) is specific activity calculated according to formula (2) at the time t_(q).

Equations of type (1) are formulated for each measured isotope of krypton and xenon; in the result, SLAE is formed

$\begin{matrix} {{{\begin{bmatrix} {a_{11}\left( {t_{q},t} \right)} & {a_{12}\left( {t_{q},t} \right)} & \ldots & {a_{1\; m}\left( {t_{q},t} \right)} \\ {a_{21}\left( {t_{q},t} \right)} & {a_{22}\left( {t_{q},t} \right)} & \ldots & {a_{2\; m}\left( {t_{q},t} \right)} \\ \vdots & \vdots & \ddots & \vdots \\ {a_{n\; 1}\left( {t_{q},t} \right)} & {a_{n\; 2}\left( {t_{q},t} \right)} & \ldots & {a_{nm}\left( {t_{q},t} \right)} \end{bmatrix} \cdot \begin{bmatrix} {\rho \; N_{1}} \\ {\rho \; N_{2}} \\ \vdots \\ {\rho \; N_{m}} \end{bmatrix}} = \begin{bmatrix} {{\overset{\sim}{A}}_{1}(t)} \\ {{\overset{\sim}{A}}_{2}(t)} \\ \vdots \\ {{\overset{\sim}{A}}_{n}(t)} \end{bmatrix}},} & (4) \end{matrix}$

wherein unknown contributions of radioactive sources ρN_(j) to total activity of krypton and xenon isotopes are to be defined.

The first step in solving the task of identifying RNG is to define separation time t_(q) for krypton and xenon isotopes.

A time interval to which separation time belongs may be found by completing “building” relative isotope activities in different fission types “in reverse time” from the measurement time with no account taking of influence of isotopes preceding thereto and by determining intersections of lines drawn from experimental points and relative activities built with due regard to influence of preceding isotopes along the isotope decay chain.

FIG. 1 shows plots of relative activities for two xenon isotopes (Xe^(133m), Xe¹³⁵), wherein the experimental point 1 corresponds to the time of measuring activities t=12 hours after an event. In order not to complicate the Figure, only “boundary” lines are shown, which correspond to U_(f) ²³⁵ and Pu_(f) ²³⁹ (instead of 6 possible fission types: U_(f) ²³⁵, U₁₄ ²³⁵, U_(f) ²³⁸, U₁₄ ²³⁸, Pu_(f) ²³⁹ and Pu₁₄ ²³⁹).

As shown in FIG. 1, the separation time belongs to the interval from t_(H)=3 to t_(K)=4 hours after the event.

When setting a time grid within the interval [t_(H), t_(K)] and solving the system (4) for times t_(q) corresponding to the grid nodes, the time t_(q) ^(min), for which a sum of squares of discrepancies of the system (4) is minimum, is taken as the separation time.

The general algorithmic diagram enabling to find estimates {circumflex over (t)}_(q), {circumflex over (ρ)}{circumflex over (N)}_(j), â_(ij) ^(true)(t_(q),t) is shown in FIG. 2.

The second step in solving the task of identifying a nuclear explosion is to determine solution estimates {circumflex over (ρ)}{circumflex over (N)}_(j) for each fixed separation time t_(q).

At a given t_(q) the system (4) is linear relative to unknown {circumflex over (ρ)}{circumflex over (N)}_(i), j=1, 2, . . . m.

Since elements of a two-dimensional signal A={a_(ij) (t_(q),t)} _(may not) b_(e ca)l_(cu)l_(ate)d accurately (independent outputs are known with errors) and isotope activities Ã(t) are also measured with errors, we can think that elements of the two-dimensional signal A and measured activities Ã(t) are independent random quantities distributed according to the normal law with mathematical expectations equal to q_(ij) ^(true)(t_(q),t) and Ã_(i) ^(true)(t) and variances equal to σ²(a_(ij)(t_(q),t)) and σ²(Ã_(i)(t)), respectively:

a _(ij)(t _(q) ,t)=a _(ij) ^(true)(t _(q) ,t)±ε_(ij),

Ã _(i)(t)=Ã_(i) ^(true)(t)±δ_(i),   (5)

where:

a_(ij) ^(true)(t_(q),t), Ã_(i) ^(true)(t) are true values of specific and measured activities of isotopes (which are unknown to us;

ε_(ij) are errors in determining specific activities a_(ij)(t_(q),t);

δ_(i) are errors in measuring activities Ã_(i)(t) of RNG in the atmosphere.

In order to take into account errors both in measured activities Ã_(i)(t) and in elements of the two-dimensional signal {a_(ij)(t_(q),t)}, the definition of orthogonal regression [4] is used, and, due to independence of random quantities a_(ij)(t_(q),t) and Ã_(i)(t)the signal may be written as:

$\begin{matrix} {{F_{l} = {\frac{1}{2}{\sum\limits_{i = 1}^{n}\begin{pmatrix} {\frac{\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)} \cdot \left( {\rho \; N_{j}} \right)}}} \right)^{2}}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)} +} \\ {\sum\limits_{j = 1}^{m}\frac{\left( {{a_{ij}\left( {t_{q},t} \right)} - {a_{ij}^{true}\left( {t_{q},t} \right)}} \right)^{2}}{\sigma^{2}\left( {a_{ij}\left( {t_{q},t} \right)} \right)}} \end{pmatrix}}}},} & (6) \end{matrix}$

where:

(ρN_(j)), j= 1,m are contributions of radioactivity sources to total activity;

a_(ij) ^(true)(t_(q),t) are unknown accurate values of specific activities which estimates are made more precise in the process of finding ρN_(j), a_(ij)(t_(q),t) are specific activities calculated according to the formulae (2)-(3) over independent and cumulative outputs, which have errors, of isobaric chain elements of radioactive transformations;

Ã_(i)(t) are RNG measured in a sample.

The signal (6), apart from an unknown vector of contributions of radioactivity sources ρN_(j), j=1, 2, . . . , m, comprises also unknown true values of calculated activities a_(ij) ^(true)(t_(q),t), i=1, 2, . . . , n, j=1, 2, . . . , m, which estimates are to be found with the use of confluent analysis [3, 4].

The following conditions are to be met in the minimum point of the signal (6):

$\begin{matrix} {{{\frac{\partial F_{i}}{\partial\left( {\rho \; N_{j}} \right)}}_{{\rho \; N_{j}} = \hat{\rho \; N_{j}}} = 0},{j = 1},2,\ldots \mspace{14mu},{m;}} & (7) \\ {{{{\frac{\partial F_{l}}{\partial{a_{ij}^{true}\left( {t_{q},t} \right)}}}_{{a_{ij}^{true}{({t_{q},t})}} = {{\hat{a}}_{ij}^{true}{({t_{q},t})}}} = 0},{i = 1},2,\ldots \mspace{14mu},{n;}}{{j = 1},2,\ldots \mspace{14mu},m}} & (8) \end{matrix}$

A structure chart of the minimum signal (6) is shown in FIG. 3; it corresponds to Block 1 of the structure chart shown in FIG. 2.

In spite of linearity at a fixed t_(q) in the equation systems (7)-(8), the task is incorrect from the point of calculations due to poor conditionality of the system (7). A relation between the maximum and minimum proper numbers in the matrix of the system (7) can be in the order of 1026. Therefore, in order to solve it, it is necessary to apply specific techniques, and in the claimed method these are methods of multi-criterion mathematical programming, wherein it is not required to determine a value for the regularization parameter, as in other methods for solving incorrect tasks.

As the first step at â_(ij) ^(true)(t_(q),t)=(t_(q),t) the SLAE (7) is solved by methods of multi-criterion mathematical programming (a method for compressing an allowable value area, target programming) and the first approximation of the ({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})_(l) estimate is found. In order to obtain estimates of true values â_(ij) ^(true)(t_(q),t) at a given value t_(q), at each step of obtaining estimates {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))}, j=1, 2, . . . , m the condition (8) [3] is used, which results in additionally solving n systems of linear equations with m unknowns of the following type:

${{{\sum\limits_{r = 1}^{m}{\frac{\left( \hat{\rho \; N_{r}} \right)\left( \hat{\rho \; N_{v}} \right)}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)}{{\hat{a}}_{ir}^{true}\left( {t_{q},t} \right)}}} + \frac{{\hat{a}}_{iv}^{true}\left( {t_{q},t} \right)}{\sigma^{2}\left( a_{iv} \right)}} = {\frac{a_{iv}\left( {t_{q},t} \right)}{\sigma^{2}\left( {a_{iv}\left( {t_{q},t} \right)} \right)} + \frac{\left( \hat{\rho \; N_{v}} \right){{\overset{\sim}{A}}_{i}(t)}}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)}}},\mspace{20mu} {i = 1},2,\ldots \mspace{14mu},n,{v = 1},2,\ldots \mspace{14mu},{m.}$

The thus obtained value estimates should comply with the natural condition, i.e., belong to the uncertainty area of measured values a_(ij)(t_(q),t)

|a _(ij)(t _(q) ,t)−â _(ij) ^(true)(t _(q) ,t)|≦3σ(q _(ij)(t _(q) t)), i=1, 2, . . . , n, j=1, 2, . . . , m.

If this condition is not met, then â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m, which do not comply with this inequality, should be replaced by values of the nearest boundary points. Due to this, values of the signal F_(l) may be increased on new constant values, as compared with those at the previous step of the iteration process, which will result in reducing the speed of iteration process convergence or in appearing oscillations.

In order values of the functional are not increased after recalculating estimates â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m the sets of estimates â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m, for which the corresponding components of the functional F_(l) are increased in comparison to their values at the previous iteration, should be replaced by the corresponding values for the previous step.

After adjusting estimates of true values â_(ij) ^(true)(t_(q),t), a next approximation {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} to solution for ρN_(j), j=1, 2, . . . , m is found by methods of multi-criterion mathematical programming instead of the regularization by A. N. Tikhonov, as is used in the analogous solution.

The criterion for stopping the algorithm is an insignificant difference between values of the signal F_(l) and components of the vector ρN_(j) in adjacent iterations, i.e., fulfillment of the following inequalities:

$\begin{matrix} {{{\frac{\left( \hat{\rho \; N_{j}} \right)_{l} - \left( \hat{\rho \; N_{j}} \right)_{l - 1}}{\left( \hat{\rho \; N_{j}} \right)_{l}}} < \gamma_{1}},{{\frac{F_{l} - F_{l - 1}}{F_{l}}} < \gamma_{2}},} & (9) \end{matrix}$

where:

({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})_(l) is next approximation to a solution at the 1st iteration;

γ₁, γ² are some numbers (small decimal fractions, e.g., 0.001) defining the accuracy of calculating values for estimates of {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))}.

When solving by methods of multi-criterion mathematical programming:

1) a two-criterion task for mathematical programming is formed:

$\begin{matrix} {{J_{1} = {{\sum\limits_{i = 1}^{n}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)}\left( {\rho \; N_{j}} \right)}}} \right)^{2}}->\min\limits_{\rho \; N_{i}}}},{J_{2} = {{\sum\limits_{j = 1}^{m}{\rho \; N_{j}}}->\min\limits_{\rho \; N_{j}}}},} & (10) \end{matrix}$

with the limitations ρN_(j)≧0, j→ 1,m.

Here, J₁ is a signal regulating a sum of squares of discrepancies for the equation system of type (4), which ensures matching a solution estimate {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))}, j= 1,m to measured activities of krypton and xenon isotopes; J₂ is a signal forming a solution type.

2) Using the threshold optimization method or target programming, the algorithm goes from the two-criterion task of mathematical programming (10) to a single-criterion task by means of transferring all the above functionals, except for one, to the limitation conditions.

The threshold optimization method (or the method of e-limitations) results in various possible combinations of target functions and limitations. The algorithm applies their following types:

$\begin{matrix} {{{{\min\limits_{\rho \; N_{j}}{\sum\limits_{i = 1}^{n}{\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)} \cdot \left( {\rho \; N_{j}} \right)}}} \right)^{2}\pi \; p\; {\sum\limits_{j = 1}^{m}{\rho \; N_{j}}}}}} \leq \delta};}{{{\rho \; N_{j}} \geq 0},{{j = \overset{\_}{1,m}};}}} & (11) \\ {{{{\min\limits_{\rho \; N_{i}}{\sum\limits_{j = 1}^{m}{\rho \; N_{j}\pi \; p\; {\sum\limits_{i = 1}^{m}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)} \cdot \left( {\rho \; N_{j}} \right)}}} \right)^{2}}}}} \leq \beta};}{{{\rho \; N_{j}} \geq 0},{j = {\overset{\_}{1,m}.}}}} & (12) \end{matrix}$

The task (11) is for quadratic programming, the task (12) is for non-linear programming.

Estimates of the right sides of the limitations δ and β may be obtained at independent minimization of the functional J1 and J2 with the limitations ρN_(j) Z≧0, j= 1,m. Here, any mathematical programming method may be used.

Target programming includes two solution models—the Archimedean model and the model with priorities.

If the Archimedean model is used, all target functions are transferred into limitations, and a weighed sum of their deviations from the limitations is minimized:

$\begin{matrix} {{{{\min\limits_{d_{1},d_{2}}{\left\{ {- \left( {{w_{1}d_{1}} + {w_{2}d_{2}}} \right)} \right\} \mspace{14mu} {at}\mspace{14mu} {J_{1}\left( {\rho \; N_{j}} \right)}}} + d_{1}} \leq \beta},{{{J_{2}\left( {\rho \; N_{j}} \right)} + d_{2}} \leq \delta},} & (13) \end{matrix}$

where:

w_(i) are weighting coefficients,

${{\sum\limits_{i = 1}^{2}w_{i}} = 1};$

are deviations from limitations.

If the model with priorities is used, target functions are successively transferred to limitations, and deviations of target function values from limitations are minimized. Here a deviation value d_(i), as found in this step, is used as an optimal deviation in the next step i+1:

$\begin{matrix} {{{{{{step}\mspace{14mu} 1\text{:}\mspace{14mu} {\min\limits_{d_{1}}{\left( {- d_{1}} \right)\mspace{14mu} {at}\mspace{14mu} {J_{1}\left( {\rho \; N_{j}} \right)}}}} + d_{1}} \leq \beta};}{{{step}\mspace{14mu} 2\text{:}\mspace{14mu} {\min\limits_{d_{2}}{\left( {- d_{2}} \right)\mspace{14mu} {at}\mspace{14mu} {J_{1}\left( {\rho \; N_{j}} \right)}}}} + {d_{1{onm}}{_{d_{1{onm}} = d_{1}}{{\leq \beta},{{{J_{2}\left( {\rho \; N_{j}} \right)} + d_{2}} \leq {\delta.}}}}}}} & (14) \end{matrix}$

For identifying a nuclear explosion by a small number of xenon isotopes the fact is used that output of such isotopes poorly depends on neutron energies.

This is illustrated in FIG. 4 that shows plots of averaged cumulative outputs of isotopes Xe133 and Xe135, depending on proportions of cumulative outputs corresponding to fission spectrum neutrons and neutrons with energy of 14 MeV. It is known that cumulative outputs have errors up to 5%.

It can be seen in FIG. 4 that averaged output values (values corresponding to equal proportions of fission spectrum neutrons and neutrons with energy of 14 MeV are shown as examples) are within these errors.

Combinations of two U²³⁵ fission types (by fission spectrum neutrons and neutrons with energy of 14 MeV) into one fission type and of two Pu239 fission types (by fission spectrum neutrons and neutrons with energy of 14 MeV) into one fission type are used for identifying by 2-4 measured isotopes, which results in reducing the number of identifiable fission types (two instead of four).

Here, specific activity is calculated as follows:

1) for fissionable U²³⁵—according to the formula (2) with the independent output vector η^(U) ²³⁵ =c₁θ^(U) ^(f) ²³⁵ +(1=c₁)η^(U) ¹⁴ ²³⁵, where η^(u) ^(f) ²³⁵ and η^(u) ¹⁴ ²³⁵ are independent outputs of elements of isobaric chains when U²³⁵ is fissioned by fission spectrum neutrons and neutrons with energy of 14 MeV, respectively; c₁ is the parameter for taking into account a proportion of independent outputs η^(U) ^(f) ²³⁵ and η^(U) ¹⁴ ²³⁵ in their sum c₁ ∈ [0,1].

2) for fissionable Pu²³⁹—according to the formula (2) with the independent output vector η^(Pu) ²³⁹ =c₂η^(Pu) ^(f) ²³⁹+(1−c₂)η^(Pu) ¹⁴ ²³⁹ , where η^(Pu) _(f) ²³⁹ and η^(Pu) ¹⁴ ²³⁹ are independent outputs of elements of isobaric chains when Pu²³⁹ is fissioned by fission spectrum neutrons and neutrons with energy of 14 MeV, respectively; c₂ is the parameter for taking into account a proportion of independent outputs η^(Pu) ^(f) ²³⁹ and η^(Pu) ¹⁴ ^(ds 239) in their sum and ⁹ c₂ ┌ [0,1].

After setting by c₁ and c₂ a two-dimensional grid with the pitches Δc₁ and Δc₂, respectively, and finding the signal minimum (6) for different c₁ and c₂, contributions of sources ρN^(U) ²³⁵ and ρN^(Pu) ²³⁹ to the total activity of krypton and xenon isotopes, at which a sum of squares of the system (2) discrepancies is minimum, are taken as true.

Thus, the proposed method for identifying parameters of a nuclear explosion can be carried out as follows:

1. Activities Ã_(i)(t) of individual krypton and xenon isotopes (Kr^(83m), Kr^(85m), Kr⁸⁵, Kr⁸⁸, Xe^(131m), Xe^(133m), Xe¹³³, Xe¹³⁵) are measured in the atmosphere (t is measurement time), where i= 1,n, n is the number of measured isotopes.

2. Dependencies of isotope activity relations on time are built with no account taken of separation from the appearance time of an event (nuclear explosion) to the time of measuring activities Ã_(i)(t) of krypton and xenon isotopes (for example, Kr^(85m)/Xe¹³⁵) for different types of fissionable materials (U_(f) ²³⁵,U₁₄ ²³⁵,Pu_(f) ²³⁹,Pu₁₄ ²³⁹).

3. Relations of isotopes (for example, Kr^(85m)/Xe¹³⁵) are determined in “reverse time” from measured experimental points until they match relative activity values for the same isotopes with no account taken of separation.

4. A time interval [t_(H),t_(K)] is determined by matching values, in which the krypton and xenon isotopes separated from their predecessors along isobaric chains of radioactive transformations.

5. A time grid with the pitch Δt is set in the interval [t_(H),t_(K)].

6. A two-dimensional signal of specific activities of krypton and xenon isotopes is formed and memorized for each node of the grid

$\left\{ {{{a_{ij}\left( {t_{q},t} \right)} = \begin{Bmatrix} {a_{ij}(t)} & {{{{\pi \; p\; t} < t_{q}},}\;} \\ {a_{ij}^{t_{q}}\left( {t_{q},t} \right)} & {{{\pi \; p\; t} > t_{q}},} \end{Bmatrix}},{i = \overset{\_}{1,n}},{j = \overset{\_}{1,m}}} \right.$

wherein lines correspond to particular isotopes and columns correspond to particular fission types; in is a number of considered fission types. Here a_(ij)(t) corresponds to the specific activity value of i-isotope for j-type of fission, which activity is calculated until the separation time t_(q) according to the formula (2), and a_(ij) ^(t) ^(q) (t_(q),t) is calculated at the time t exceeding the separation time t_(q) according to the formula (3).

7. Considering that the measured activity values Ã_(i)(t) and the elements of the two-dimensional signal {a_(ij)(t_(q),t)} are statistically independent quantities distributed according to the normal law with mathematical expectations Ã_(i) ^(true)(t) and a_(ij) ^(true)(t_(q),t) and mean square deviations σ(Ã_(i)(t)) and σ(a_(ij)(t_(q),t)), respectively, a one-dimensional signal F_(l) is formed according to the formula (6).

8. Numbers γ₁, γ₂, which characterize accuracy of estimating {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} and â_(ij) ^(true)(t_(q),t) are set.

9. A two-criterion task for mathematical programming

${J_{1} = {{\sum\limits_{i = 1}^{n}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{a_{ij}^{true}\left( {t_{q},t} \right)} \cdot \left( {\rho \; N_{j}} \right)}}} \right)^{2}}->\min\limits_{\rho \; N_{j}}}},{J_{2} = {{\overset{m}{\sum\limits_{j = 1}}{\rho \; N_{j}}}->\min\limits_{\rho \; N_{j}}}},$

is formed in the first iteration, taking a_(ij) ^(true)(t_(q),t)=a_(ij)(t_(q),t) (and in the following iterations a_(ij) ^(true)(t_(q),t)=â_(ij) ^(true)(t_(q),t)) with the limitations ρN_(j)≧0, j= 1,m.

10. Using the method for compressing allowable value area or target programming, the algorithm goes from the mathematical programming two-criterion task according to the formulae (11)-(14) to a one-criterion task by means of transferring all the above functionals, except for one, to the limitation conditions.

11. Using methods of quadratic programming, non-linear programming, target programming (the Archimedean model and the model with priorities), the first approximation ({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})_(l) of estimated activity contributions of individual sources to the total activity.

12. After obtaining the first approximation ({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})_(l), the elements of the two-dimensional signal {a_(ij)(t_(q),t)} are specified more precisely. For this, n systems of linear equations with in unknowns are solved:

${{{\sum\limits_{r = 1}^{m}{\frac{\left( {\rho \; N_{r}} \right)\left( {\rho \; N_{r}} \right)}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)}{{\hat{a}}_{ir}^{true}\left( {t_{q},t} \right)}}} + \frac{{\hat{a}}_{iv}^{true}\left( {t_{q},t} \right)}{\sigma^{2}\left( {a_{iv}\left( {t_{q},t} \right)} \right)}} = {\frac{a_{iv}\left( {t_{q},t} \right)}{\sigma^{2}\left( {a_{iv}\left( {t_{q},t} \right)} \right)} + \frac{\left( {\rho \; N_{v}} \right){{\overset{\sim}{A}}_{i}(t)}}{\sigma^{2}\left( {{\overset{\sim}{A}}_{i}(t)} \right)}}},\mspace{20mu} {i = 1},2,\ldots \mspace{14mu},n,{v = 1},2,\ldots \mspace{14mu},m,$

where the estimate approximation ({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})^(l) found in the first iteration is substituted for (ρN_(j)).

13. Then the algorithm checks whether the new values â_(ij) ^(true)(t_(q),t) satisfy the (natural area of uncertainty of the elements a_(ij)(t_(q),t):

|a _(ij)(t _(j) ,t)−â _(ij) ^(true)(t _(q) ,t)|≦3σ(a _(ij)(t _(q) ,t)), i=1, 2, . . . , n, j=1, 2, . . . , m.

If this condition is not met, then â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m, which do not satisfy this inequality, are replaced by values of the nearest boundary points.

Due to this, values of the signal F_(l) may be greater on new values of variables â_(ij) ^(true)(t_(q),t) in comparison to the previous step of the iteration process, which results in reducing the speed of convergence of the iteration process and appearing oscillations.

In order values of the signal F_(l) are not increased after recalculating estimates â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m, the sets of estimates â_(ij) ^(true)(t_(q),t), j=1, 2, . . . , m, for which the corresponding components of the functional F_(l) are increased in comparison to their values at the previous iteration, should be replaced by the corresponding values for the previous step.

14. The operations mentioned in Items 9-13 should be repeated until the following conditions are met:

${{\frac{\left( \hat{\rho \; N_{j}} \right)_{l} - \left( \hat{\rho \; N_{j}} \right)_{l - 1}}{\left( \hat{\rho \; N_{j}} \right)_{l}}} < \gamma_{1}},{{\frac{F_{l} - F_{l - 1}}{F_{l}}} < \gamma_{2}},$

15. For identifying a nuclear explosion by a small number of isotopes (2-4 isotopes) a two-dimensional grid is set over c₁ and c₂ with the pitches Δc₁ and Δc₂, respectively, where c₁ and c₂ are weights for summing up the independent and cumulative outputs of the RTIC elements.

16. Vectors of independent outputs η^(U) ²³⁵ =c₁η^(U) ^(f) ²³⁵ +(1−c₁)η^(U) ¹⁴ ²³⁵ and η^(Pu) ²³⁹ =c₂η^(Pu) ^(f) ²³⁹+(1−c₂)η^(Pu) ¹⁴ ²³⁹ are calculated for each value of c₁ and c₂.

17. A two-dimensional signal of specific activity calculated and memorized for each pair of the vectors η^(U) ²³⁵ and η^(Pu) ²³⁹ .

18. The operations mentioned in Items 6-13 are performed for each two-dimensional signal, until the condition from Item 14 is met.

19. Among the signal estimates ({circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))})_(c) ¹ _(c) ² , which are found after solving the task of identifying a nuclear explosion, at the corresponding two-dimensional signal {a_(ij)(t_(q),t)}_(c) ¹ _(ds 2), that estimate is selected at which a value of

$\sum\limits_{i = 1}^{n}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{\left( {a_{ij}^{true}\left( {t_{q},t} \right)} \right)_{c_{1}c_{2}}\left( \hat{\rho \; N_{j}} \right)_{c_{1}c_{2}}}}} \right)^{2}$

is minimal.

20. As the separation time tq is taken at which the relation

$\sum\limits_{i = 1}^{n}{\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{{{\hat{a}}_{ij}^{true}\left( {t_{q},t} \right)}\left( \hat{\rho \; N_{j}} \right)}}} \right)^{2}/{\sum\limits_{j = 1}^{m}\left( \hat{\rho \; N_{j}} \right)^{2}}}$

is minimal.

Shown in the general block diagram of the algorithm (FIG. 2), Block “Introducing separation interval [t_(H),t_(K)] and pitch Δt within this interval” corresponds to Items 4, 5; Block “Minimizing signal f_(l)(t_(q),{circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))},â_(ij) ^(true)) by t_(q) corresponds to Item 20.

The main part of the algorithm is put into Block 1 “Searching estimates for {circumflex over (ρ)}{circumflex over (n)}{circumflex over (ρn_(j))},â_(ij) ^(true)” which diagram is shown in FIG. 3.

Block “Input of initial data t_(q),Ã_(i)(t),σ(Ã_(i)(t)),σ(a_(ij)(t_(ij),t)),γ₁,γ₂ corresponds to Items 1, 7, 8 of the method;)

Block “Calculation of elements of two-dimensional signal {a_(ij)(t_(q),t)} for case of instant fission” corresponds to Items 6, 17 of the method;

Block “Grid by c₁ and c₂. Forming averaged independent and cumulative outputs” corresponds to Items 15, 16 of the method;

Block “Initializing â_(ij) ^(true)(t_(q),t)=a_(ij)(t_(q),t)” corresponds to Item 9 of the method. The operations mentioned in Items 9-11 of the method are performed in Block “Searching for {circumflex over (ρ)}{circumflex over (N)}{circumflex over (ρN_(j))} approximation by methods of multi-criterion mathematical programming”;

Block “Solving systems of linear equations for obtaining next â_(ij) ^(true)(t_(q),t) approximation realizes Item 12, and Blocks “|â_(ij) ^(true)(t_(q),t)−a_(ij)(t_(q),t)|≦3σ(a_(ij)(t_(q),t))”, “Replacing estimates â_(ij) ^(true)(t_(q),t) by boundary point values”, “Replacing estimates â_(ij) ^(true)(t_(q),t) by values of previous iteration” and the conditions binding them correspond to Item 13 of the method. lock “Output condition met” realizes Item 14 of the method, and Block “Selecting optimal solution for averaged outputs” realizes Item 19.

INDUSTRIAL APPLICABILITY OF THE PROPOSED METHOD

Implementation of the proposed method was simulated on a personal computer having 2.4 GHz Intel Celeron processor, RAM 768 MB and the Matlab 7.0 mathematical package.

The simulation was directed to a situation of sample taking 6 days after an explosion and to measuring activities of 5 isotopes (Kr^(85m), Xe^(131m), Xe^(133m), Xe¹³³, Xe¹³⁵). Results were calculated on the condition that separation time is supposed to be known and equal to 3 hours after the event.

Values of measured activities were additively “made noisy” with Gaussian noise with root-mean-square deviation equal to 5% of their “true” value.

The following combinations of fission types (possible sets of ρN_(j) variables):

1) U_(th) ²³⁵+Xe¹³³ background (two unknown sources);

2) U_(th) ²³⁵+U_(f) ²³⁵+U₁₄ ²³⁵ (three unknown sources);

3) U_(th) ²³⁵+Pu_(f) ²³⁹+Pu₁₄ ²³⁹ (three unknown sources);

4) U_(f) ²³⁵+U₁₄ ²³⁵ (two unknown sources);

5) U_(f) ²³⁵+U₁₄ ²³⁵ +Xe¹³³ background (three unknown sources);

6) Pu_(f) ²³⁹+Pu₁₄ ²³⁹ (two unknown sources);

7) Pu_(f) ²³⁹+Pu₁₄ ²³⁹+Xe¹³³ background (three unknown sources);

8) U_(f) ²³⁵+U₁₄ ²³⁵+Pu_(f) ²³⁹+Pu₁₄ ²³⁹ (four unknown sources);

9) U_(f) ²³⁵+U₁₄ ²³⁵+Pu_(f) ²³⁹+Pu₁₄ ²³⁹+Xe¹³³ background (five unknown sources), where U_(th) ²³⁵ is reactor emission (data on reactors is taken from literature).

The true solution is the combination 4, the relative contribution of U_(f) ²³⁵ is 100, the relative contribution of U₁₄ ²³⁵ is 100. The simulation results are shown in Table 1.

The line “Method for solving” states 4 methods for solving the task of identifying a nuclear explosion (quadratic programming, non-linear programming, Archimedean model, model with priorities), as proposed in this method, which were compared to the solving method that is used in the analogous algorithm (Tiknonov's regularization).

The line “Number of fission type combination” indicates the number of fission type combination that ensures the least sum of squares of discrepancies in the equation system (4) for corresponding solving methods (Tikhonov's regularization, methods of multi-criterion mathematical programming) from the 9 combinations in total.

The line “Order of conditionality number for system matrix” indicates orders of conditionality numbers for the matrix of the system (4), as corresponds to the fission type combination indicated in the second line.

The line “Solution estimate” indicates estimates for fission type contributions present in fission type combinations indicated in the second line. For example, the combination 2 is the best for the Tikhonov's regularization from the point of the sum of discrepancy squares. Three fission types correspond to this combination: U_(th) ²³⁵, which calculated relative contribution is 35.99; U_(f) ²³⁵, which calculated relative contribution is 42.97, and U₁₄ ²³⁵, which calculated relative contribution is 110.27. The same applies to the other solving methods also.

The line “Sum of discrepancy squares” indicates values of discrepancy square sum for the equation system (4), which are calculated for fission type combinations and their estimated contributions to the total activity of the krypton and xenon isotopes.

The line “Time of using algorithm, min” indicates time in minutes required for obtaining contribution estimates by the respective method.

TABLE 1 Results of solving the identification task by various methods (radioactivity source are U₂₃₅ ^(f) and U₂₃₅ ¹⁴, true solution is 100 and 100) Tikhonov's Quadratic Non-linear Archimedean Model with Solving method regularization programming programming model priorities Number of 2 9 9 9 9 fission type combination Order of 10⁶ 10¹⁴ 10¹⁴ 10¹⁴ 10¹⁴ conditionally number for system matrix Solution 35.99 84.94 43.53 84.94 84.94 estimate 42.97 108.18 120.76 108.18 108.18 110.27 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 10.99 0.00 10.98 10.99 0.00 0.00 0.00 0.00 Sum of 751.39 74.97 5096.43 74.97 74.97 discrepancy squares Time of using 1.46 9.28 11.95 13.40 19.78 algorithm, min

It can be seen in Table 1 that the Tikhonov's regularization method gives a negative result, the solution comprises significant relative contribution of a nuclear reactor (the true solution lacks it).

A positive solution is obtained by the methods of multi-criterion programming (quadratic, non-linear, target programming (Archimedean model and model with priorities)), when the additional condition of non-negativity of variables is used.

The optimal result corresponds to the ninth combination of fission types. This does not contradict to the true solution, since contributions of those fission types that are not present in the true solution are insignificant (most of them are zero).

Thus, the advantages of the proposed method are:

Increased efficiency and accuracy of identification of nuclear explosion parameters due to: simultaneously checking various fission type combinations by methods of multi-criterion mathematical programming, which do not require calculation of the regularization parameter, using additional limitations applied to a solution and corresponding the task physical statement (non-negativity, boundedness of solution), using different target functions simultaneously and combining two fission types of one material by neutrons from two energy groups into a single fission type, i.e., reducing the number of unknowns in equation systems. 

1. A method for identifying a nuclear explosion by radioactive isotopes of krypton and xenon, characterized by measuring, at time t after an event, signals Ã_(i)(t) describing a change in total activity of each isotope in the atmosphere near a measuring station; building relations between isotope activities and time with no account taken of separation and relations of isotope activities drawn from measured points in reverse time for all fission types under consideration; determining a separation interval [t_(H), t_(K)]; setting a time grid within the separation interval [t_(H), t_(K)]; forming fission type combinations; calculating t_(q) for each grid node within the separation interval and for each combination a two-dimensional signal A describing values of “specific” activities of each isotope against the separation time t_(q) and measuring time t and potential source (of fissionable material and neutron energy); setting root-mean-square values σ(Ã_(i)(t)) for errors in measured signals Ã_(i)(t) and errors σ(a_(ij)(t_(q),t)) in the two-dimensional signal elements {a_(i:j)(t_(q),t); setting quantities γ₁,γ₂ defining accuracy of signal estimates {circumflex over (p)}{circumflex over (N)}_(j); identifying a nuclear explosion by means of forming a signal J₁ defining accuracy of a solution, which is obtained from a sum of squares of difference between signals Â_(i)(t) and products obtained by multiplying lines of a two-dimensional signal {a_(ij)(t_(q),t)} by values of {circumflex over (p)}{circumflex over (N)}î{circumflex over (pNi_(j))} signals, and a signal J₂ defining a type of a {circumflex over (p)}{circumflex over (N)}{circumflex over (pN_(j))} signal; forming limitation signals and a target signal (target function) from the signals J₁ and J₂ in the set combinations; and finding, with the use of an iteration process and subject to specifying at each iteration, two-dimensional signal values {a_(ij)(t_(q),t)} for elements of the {circumflex over (p)}{circumflex over (p_(j))} signal; determining, according to the elements of the {circumflex over (p)}{circumflex over (N)}{circumflex over (pN_(j))} signal, point estimates foractivity contributions of each possible fission type to total activity; selecting an optimal combination of fission types according to a sum of discrepancy squares.
 2. A method according to claim 1, characterized by combining two types of fission of U²³⁵ and Pu²³⁹ by neutrons of different energy groups into one fission type η^(U) ²³⁵ and η^(Pu) ²³⁹ , respectively, by summing up with weights of independent and cumulative outputs of isobaric chain elements corresponding to different fission types; forming a two-dimensional grid by weights c₁ and c₂; calculating, for each pair of weights (c₁, c₂), elements of a two-dimensional signal of specific activity {a_(ij)(t_(q),t)}_(c) ₁ _(c) ² ; calculating, according to claim 1, estimates for signals {circumflex over (p)}{circumflex over (N)}j)_(c) ₁ _(c) ₂ ; selecting a ({circumflex over (p)}{circumflex over (N)}j)_(c) ₁ _(c) ₂ estimate for which a value of $\sum\limits_{i = 1}^{n}\left( {{{\overset{\sim}{A}}_{i}(t)} - {\sum\limits_{j = 1}^{m}{\left( {a_{ij}^{true}\left( {t_{q},t} \right)} \right)_{c_{1}c_{2}}\left( \hat{\rho \; N_{j}} \right)_{c_{1}c_{2\;}}}}} \right)^{2}$ is minimal. 